## This R file Runs the Bayesian Factor Analysis presented in paper#
## Code written by Jonathan Slapin#
## Final version: October 12, 2011#
#
rm(list=ls())#
#
library(calibrate)#
library(foreign)#
library(MCMCpack)#
#
# Set working directory#
setwd("/Users/jslapin/Dropbox/RTA Exit Options/rio/replication files")#
#setwd("/Users/jgray/Documents/My Dropbox/RTA Exit Options/rio/replication files")#
#
# Read in final dataset from Stata#
rta<-as.data.frame(read.dta("rta_meanXsal.dta"))#
#
# Subset data to eliminate highly missing RTAs#
miss<-ifelse(is.na(rta[,2:ncol(rta)])==T,1,0)#
colmiss<-c(0,apply(miss,2,sum))#
rta<-rta[,colmiss<20]#
rowmiss<-apply(miss,1,sum)#
rta<-rta[rowmiss<20,]#
rownames(rta)<-rta[,1]#
rta<-rta[,-1]#
#
# Bayesian factor analysis#
post1dimshort<- MCMCfactanal(~Ability + Ambition + Delegation + DSM + Enforcement + Goods + IntraPTA + Legalization + Monitoring +  NTBs  + Scope + Secretariat + Services +  Tariffinternal + Tradediversion,factors=1,lambda.constraints=list(Ability=list(1,"+"),Secretariat=list(1,"+"), IntraPTA=list(1,"+") ), verbose=0, store.scores=TRUE, a0=1, b0=0.15, data=rta, burnin=50000,mcmc=200000,thin=200)
rm(list=ls())#
#
library(calibrate)#
library(foreign)#
library(MCMCpack)#
#
# Set working directory#
setwd("/Users/jslapin/Dropbox/RTA Exit Options/rio/replication files")#
#setwd("/Users/jgray/Documents/My Dropbox/RTA Exit Options/rio/replication files")#
#
# Read in final dataset from Stata#
rta<-as.data.frame(read.dta("rta_meanXsal.dta"))#
#
# Subset data to eliminate highly missing RTAs#
miss<-ifelse(is.na(rta[,2:ncol(rta)])==T,1,0)#
colmiss<-c(0,apply(miss,2,sum))#
rta<-rta[,colmiss<20]#
rowmiss<-apply(miss,1,sum)#
rta<-rta[rowmiss<20,]#
rownames(rta)<-rta[,1]#
rta<-rta[,-1]
post1dimshort<- MCMCfactanal(~Ability + Ambition + Delegation + DSM + Enforcement + Goods + IntraPTA + Legalization + Monitoring +  NTBs  + Scope + Secretariat + Services +  Tariffinternal + Tradediversion,factors=1,lambda.constraints=list(Ability=list(1,"+"),Secretariat=list(1,"+"), IntraPTA=list(1,"+") ), verbose=0, store.scores=TRUE, a0=1, b0=0.15, data=rta, burnin=50000,mcmc=200000,thin=200)
names(rta)
#
rm(list=ls())#
#
library(calibrate)#
library(foreign)#
library(MCMCpack)#
#
# Set working directory#
setwd("/Users/jslapin/Dropbox/RTA Exit Options/rio/replication files")#
#setwd("/Users/jgray/Documents/My Dropbox/RTA Exit Options/rio/replication files")#
#
# Read in final dataset from Stata#
rta<-as.data.frame(read.dta("rta_meanXsal.dta"))#
#
# Subset data to eliminate highly missing RTAs#
miss<-ifelse(is.na(rta[,2:ncol(rta)])==T,1,0)#
colmiss<-c(0,apply(miss,2,sum))#
rta<-rta[,colmiss<20]#
rowmiss<-apply(miss,1,sum)#
rta<-rta[rowmiss<20,]#
rownames(rta)<-rta[,1]#
rta<-rta[,-1]
# Bayesian factor analysis#
post1dimshort<- MCMCfactanal(~Ability + Ambition + Delegation + DSM + Enforcement + Goods + IntraPTA + Legalization + Monitoring +  NTBs  + Scope + Secretariat + Services +  Tariffinternal + Tradediversion,factors=1,lambda.constraints=list(Ability=list(1,"+"),Secretariat=list(1,"+"), IntraPTA=list(1,"+") ), verbose=0, store.scores=TRUE, a0=1, b0=0.15, data=rta, burnin=50000,mcmc=200000,thin=200)
sum1dim<-summary(post1dimshort)
sum1dim
# Extract RTA position scores (phi) from factor analysis#
phimeans<-sum1dim$statistics[31:nrow(sum1dim$statistics),1]#
#
# Create Bayesian Credible Intervals#
ci<-HPDinterval(post1dimshort,prob=0.95)#
philci<-ci[31:nrow(ci),1]#
phihci<-ci[31:nrow(ci),2]
philci
phihci
ndf2<-read.dta("comparisons.dta")#
ndf2$phimeans<-bayesciphi$phimeans#
fit<-lm(phimeans~integration,data=ndf2)#
fit2<-lm(phimeans~integration,data=ndf2[-10,])#
cor(ndf2$integration,ndf2$phimeans)
ndf2<-read.dta("comparisons.dta")#
ndf2$phimeans<-bayesciphi$phimeans#
fit<-lm(phimeans~integration,data=ndf2)#
fit2<-lm(phimeans~integration,data=ndf2[-10,])#
cor(ndf2$integration,ndf2$phimeans)
ndf2<-read.dta("comparisons.dta")
ndf2$phimeans<-bayesciphi$phimeans
# Create dataframe with RTA names and Bayesian Credible Intervals for plotting#
rtanames<-cbind("ALADI","APEC","Andean Community","Asean","Bangkok Agreement","CARICOM","CEFTA","ECOWAS","EFTA","EU","GCC","Mercosur","NAFTA","OECS","SACU","SADC","SPARTECA","UEMOA")#
bayesciphi<-as.data.frame(cbind(philci,phimeans,phihci))#
rownames(bayesciphi)<-rtanames#
bayesciphi<-bayesciphi[order(phimeans,decreasing=TRUE),]
ndf2<-read.dta("comparisons.dta")#
ndf2$phimeans<-bayesciphi$phimeans#
fit<-lm(phimeans~integration,data=ndf2)#
fit2<-lm(phimeans~integration,data=ndf2[-10,])#
cor(ndf2$integration,ndf2$phimeans)
plot(ndf2$integration,ndf2$phimeans,pch=20,xlim=c(0,4.8),ylim=c(-.55,1.1),xlab="Bearce-Omori Integration Score",ylab="Expert Survey Principal Dimension")#
textxy(ndf2$integration,ndf2$phimeans,ndf2$rtaname,cx=.75)#
abline(fit)#
abline(fit2,lty=2)
plot(ndf3$oversightfunctions,ndf3$phimeans,pch=20,xlim=c(0,4.8),ylim=c(-.55,1.1), xlab="Grigorescu Transparency Score",ylab="Expert Survey Principal Dimension")#
textxy(ndf3$oversightfunctions,ndf3$phimeans,ndf3$rtaname,cx=0.75)#
abline(fit)
fit<-lm(phimeans~oversightfunctions,data=ndf2)#
fit2<-lm(phimeans~ oversightfunctions,data=ndf2[-10,])#
cor(ndf2$oversightfunctions,ndf2$phimeans)#
ndf3<-ndf2[is.na(ndf2$oversightfunctions)!=T,]
plot(ndf3$oversightfunctions,ndf3$phimeans,pch=20,xlim=c(0,4.8),ylim=c(-.55,1.1), xlab="Grigorescu Transparency Score",ylab="Expert Survey Principal Dimension")#
textxy(ndf3$oversightfunctions,ndf3$phimeans,ndf3$rtaname,cx=0.75)#
abline(fit)
fit<-lm(ndf$phimeans~ndf$V1)#
cor(ndf$V1,ndf$phimeans)#
plot(ndf$V1,ndf$phimeans,pch=20,xlim=c(0,4.8),ylim=c(-.55,1.1),xlab="McCall Smith Legalization Score",ylab="Expert Survey Principal Dimension")#
textxy(ndf$V1,ndf$phimeans,ndf$Row.names,cx=.75)#
abline(fit)
mccsmithrta<-rbind("Andean Community","CARICOM","CEFTA","ECOWAS","EFTA","EU","GCC","Mercosur","NAFTA","OECS","SACU")#
mccsmithscore<-t(as.matrix(cbind(4,1,0,3,0,4,1,2,2,2,0)))#
rownames(mccsmithscore)<-mccsmithrta#
ndf<-merge(mccsmithscore,bayesciphi, by.x=0,by.y=0)#
ndf$Row.names[1]<-"Andean Com"
fit<-lm(ndf$phimeans~ndf$V1)#
cor(ndf$V1,ndf$phimeans)#
plot(ndf$V1,ndf$phimeans,pch=20,xlim=c(0,4.8),ylim=c(-.55,1.1),xlab="McCall Smith Legalization Score",ylab="Expert Survey Principal Dimension")#
textxy(ndf$V1,ndf$phimeans,ndf$Row.names,cx=.75)#
abline(fit)
